#ifndef _MATHUTILS_H
#define _MATHUTILS_H

#include <math.h>

typedef unsigned short float16;

// This ain't a 1.5.10 float16, no sire, it's a stupid hack format where we chop 16 bits off a float.
// This choice is subject to change.
inline float16 FloatToFloat16(float x) {
	int ix;
	memcpy(&ix, &x, sizeof(float));
	return ix >> 16;
}

inline float Float16ToFloat(float16 ix) {
	float x;
	memcpy(&x, &ix, sizeof(float));
	return x;
}


#define PI 3.141592653589793f

// The stuff in this file is from all over the web, esp. dspmusic.org. I think it's all public domain.
// In any case, very little of it is used anywhere at the moment.

inline float sine(float t,float f,float ph,float fm) {
	return sinf((t*f+ph)*2*PI + 0.5f*PI*fm*(1 - sqrt(f*2)));
}

//fb := feedback (0 to 1) (1 max saw)

inline float saw(float t,float f,float ph, float fm, float fb = 1.0f)
{
	return sine(t,f,ph,fb*sine(t-1.0f,f,ph,fm));
}

//	pm := pulse mod (0 to 1) (1 max pulse)
//	pw := pulse width (0 to 1) (1 square)
inline float pulse(float t,float f,float ph,float fm,float fb,float pm,float pw) {
	return saw(t,f,ph,fm,fb) - saw(t,f,ph+0.5f*pw,fm,fb) * pm;
}

// Calculate pseudo-random 32 bit number based on linear congruential method. 
void SetSeed(unsigned int seed);
unsigned int GenerateRandomNumber();
inline float GenerateRandomFloat01() {
	return (float)((double)GenerateRandomNumber() / 0xFFFFFFFF);
}
inline float GenerateRandomSignedFloat() {
	return (float)((double)GenerateRandomNumber() / 0x80000000) - 1.0f;
}


inline float GaussRand()
{
	float R1 = GenerateRandomFloat01();
	float R2 = GenerateRandomFloat01();

	float X = sqrtf( -2.0f * logf(R1)) * cosf(2.0f * PI * R2);
	if (X > 4.0f) X = 4.0f;
	if (X < -4.0f) X = -4.0f;
	return X;
}

// Accuracy unknown
inline double atan_fast(double x) {
    return (x / (1.0 + 0.28 * (x * x)));
}




/*
#define PINK_NOISE_NUM_STAGES 3

class PinkNoise {
public:
	PinkNoise() {
		srand ( 0xcdebae); // initialize random generator
		clear();
	}

	void clear() {
		for( size_t i=0; i< PINK_NOISE_NUM_STAGES; i++ )
			state[ i ] = 0.0;
	}

	float tick() {
		static const float RMI2 = 2.0 / float(RAND_MAX); // + 1.0; // change for range [0,1)
		static const float offset = A[0] + A[1] + A[2];

		// unrolled loop
		float temp = float( rand() );
		state[0] = P[0] * (state[0] - temp) + temp;
		temp = float( rand() );
		state[1] = P[1] * (state[1] - temp) + temp;
		temp = float( rand() );        
		state[2] = P[2] * (state[2] - temp) + temp;
		return ( A[0]*state[0] + A[1]*state[1] + A[2]*state[2] )*RMI2 - offset;
	}

protected:
	float state[ PINK_NOISE_NUM_STAGES ];
	static const float A[ PINK_NOISE_NUM_STAGES ];
	static const float P[ PINK_NOISE_NUM_STAGES ];
};

const float PinkNoise::A[] = { 0.02109238, 0.07113478, 0.68873558 }; // rescaled by (1+P)/(1-P)
const float PinkNoise::P[] = { 0.3190,  0.7756,  0.9613  };

double interpolate(double* wavetable, int nMask, double location)
{
	// 5-point spline

	int nearest_sample = (int) location;
	double x = location - (double) nearest_sample;

	double p0=wavetable[(nearest_sample-2)&nMask];
	double p1=wavetable[(nearest_sample-1)&nMask];
	double p2=wavetable[nearest_sample];
	double p3=wavetable[(nearest_sample+1)&nMask];
	double p4=wavetable[(nearest_sample+2)&nMask];
	double p5=wavetable[(nearest_sample+3)&nMask];

	return p2 + 0.04166666666*x*((p3-p1)*16.0+(p0-p4)*2.0
	+ x *((p3+p1)*16.0-p0-p2*30.0- p4
	+ x *(p3*66.0-p2*70.0-p4*33.0+p1*39.0+ p5*7.0- p0*9.0
	+ x *( p2*126.0-p3*124.0+p4*61.0-p1*64.0- p5*12.0+p0*13.0
	+ x *((p3-p2)*50.0+(p1-p4)*25.0+(p5-p0)*5.0)))));
}
*/
/*
t := time (0 to ...)
ph := phase shift (0 to 1)
fm := freq mod (0 to 1)
*/

inline float ThirdInterp(const float x, const float L1, const float L0, const
						 float H0, const float H1) {
	return
		L0 +
		.5f*
		x*(H0-L1 +
		x*(H0 + L0*(-2) + L1 +
		x*( (H0 - L0)*9 + (L1 - H1)*3 +
		x*((L0 - H0)*15 + (H1 -  L1)*5 +
		x*((H0 - L0)*6 + (L1 - H1)*2 )))));
}	

// linear -> dB conversion
inline float lin2dB(float lin) {
	const float LOG_2_DB = 8.6858896380650365530225783783321f;	// 20 / ln( 10 )
	return log(lin) * LOG_2_DB;
}

// dB -> linear conversion
inline float dB2lin(float dB) {
	const float DB_2_LOG = 0.11512925464970228420089957273422f;	// ln( 10 ) / 20
	return exp(dB * DB_2_LOG);
}

#endif
